需要自行下载安装一些必要的R包,主要是scRNAseq,然后是其它一些辅助包用来探索这个数据集。
因为大量学员在中国大陆,通常不建议大家使用下面的R包安装方法,建议是切换镜像后再下载R包。参考:http://www.bio-info-trainee.com/3727.html
if (!requireNamespace("Rtsne"))
install.packages("Rtsne")
if (!requireNamespace("FactoMineR"))
install.packages("FactoMineR")
if (!requireNamespace("factoextra"))
install.packages("factoextra")
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
if (!requireNamespace("scater"))
BiocManager::install("scater")
if (!requireNamespace("scRNAseq"))
BiocManager::install("scRNAseq")
if (!requireNamespace("M3Drop"))
BiocManager::install("M3Drop")
if (!requireNamespace("ROCR"))
BiocManager::install("ROCR")
加载R包,前提是已经成功安装了。
rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))
library(ggplot2)
library(tidyr)
library(cowplot)
library("FactoMineR")
library("factoextra")
library("ROCR")
这个包内置的是 Pollen et al. 2014 数据集,人类单细胞细胞,分成4类,分别是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ,还有 “GW16” and “GW21” ,“GW21+3” 这种孕期细胞,理解这些需要一定的生物学背景知识,如果不感兴趣,可以略过。
这个R包大小是50.6 MB,下载需要一点点时间,先安装加载它们。
这个数据集很出名,截止2019年1月已经有近400的引用了,后面的人开发R包算法都会在其上面做测试,比如 SinQC 这篇文章就提到:We applied SinQC to a highly heterogeneous scRNA-seq dataset containing 301 cells (mixture of 11 different cell types) (Pollen et al., 2014).
不过本例子只使用了数据集的4种细胞类型而已,因为 scRNAseq 这个R包就提供了这些,完整的数据是 23730 features, 301 samples, 地址为https://hemberg-lab.github.io/scRNA.seq.datasets/human/tissues/ , 这个网站非常值得推荐,简直是一个宝藏。
这里面的表达矩阵是由 RSEM (Li and Dewey 2011) 软件根据 hg38 RefSeq transcriptome 得到的,总是130个文库,每个细胞测了两次,测序深度不一样。
library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
## SRR1275356 SRR1274090 SRR1275251 SRR1275287
## A1BG 0 0 0 0
## A1BG-AS1 0 0 0 0
## A1CF 0 0 0 0
## A2M 0 0 0 33
sample_ann <- as.data.frame(colData(fluidigm))
DT::datatable(sample_ann)
前面说到,这个数据集是130个文库,每个细胞测了两次,测序深度不一样,这130个细胞,分成4类,分别是: pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ,还有 “GW16” and “GW21” ,“GW21+3” 这种孕期细胞。
批量,粗略的看一看各个细胞的一些统计学指标的分布情况
box <- lapply(colnames(sample_ann[,1:19]),function(i) {
dat <- sample_ann[,i,drop=F]
dat$sample=rownames(dat)
## 画boxplot
ggplot(dat, aes('all cells', get(i))) +
geom_boxplot() +
xlab(NULL)+ylab(i)
})
plot_grid(plotlist=box, ncol=5 )
# ggsave(file="stat_all_cells.pdf")
很明显,他们可以有根据来进行分组,这里不再演示。 不过通常的文章并不会考虑如此多的细节,这里重点是批量,代码技巧非常值得你们学校。
因为进行了简单探索,对表型数据就有了把握,接下来可以进行一定程度的过滤,因为细节太多,这里重点是批量,代码技巧非常值得你们学校。
pa <- colnames(sample_ann[,c(1:9,11:16,18,19)])
tf <- lapply(pa,function(i) {
# i=pa[1]
dat <- sample_ann[,i]
dat <- abs(log10(dat))
fivenum(dat)
(up <- mean(dat)+2*sd(dat))
(down <- mean(dat)- 2*sd(dat) )
valid <- ifelse(dat > down & dat < up, 1,0 )
})
tf <- do.call(cbind,tf)
choosed_cells <- apply(tf,1,function(x) all(x==1))
table(sample_ann$Biological_Condition)
##
## GW16 GW21 GW21+3 NPC
## 52 16 32 30
sample_ann=sample_ann[choosed_cells,]
table(sample_ann$Biological_Condition)
##
## GW16 GW21 GW21+3 NPC
## 36 11 23 29
ct <- ct[,choosed_cells]
ct[1:4,1:4]
## SRR1274090 SRR1275287 SRR1275364 SRR1275269
## A1BG 0 0 0 0
## A1BG-AS1 0 0 0 0
## A1CF 0 0 0 0
## A2M 0 33 0 51
counts <- ct
fivenum(apply(counts,1,function(x) sum(x>0) ))
## A1CF OR8G1 LINC01003 MRPS36 YWHAZ
## 0 0 4 26 99
boxplot(apply(counts,1,function(x) sum(x>0) ))
fivenum(apply(counts,2,function(x) sum(x>0) ))
## SRR1275365 SRR1275345 SRR1275248 SRR1275273 SRR1274125
## 1566.0 3043.0 4002.0 5706.5 8024.0
hist(apply(counts,2,function(x) sum(x>0) ))
choosed_genes=apply(counts,1,function(x) sum(x>0) )>0
table(choosed_genes)
## choosed_genes
## FALSE TRUE
## 9496 16759
counts <- counts[choosed_genes,]
下面的计算,都是基于log后的表达矩阵。
dat <- log2(edgeR::cpm(counts) + 1)
dat[1:4, 1:4]
## SRR1274090 SRR1275287 SRR1275364 SRR1275269
## A1BG 0 0.000000 0 0.000000
## A1BG-AS1 0 0.000000 0 0.000000
## A2M 0 4.216768 0 3.552694
## A2M-AS1 0 0.000000 0 0.000000
dat_back <- dat
先备份这个表达矩阵,后面的分析都用得上
exprSet <- dat_back
colnames(exprSet)
## [1] "SRR1274090" "SRR1275287" "SRR1275364" "SRR1275269" "SRR1275263"
## [6] "SRR1274117" "SRR1274089" "SRR1275248" "SRR1275345" "SRR1274125"
## [11] "SRR1275300" "SRR1275294" "SRR1274122" "SRR1275277" "SRR1275293"
## [16] "SRR1274128" "SRR1275348" "SRR1275245" "SRR1274084" "SRR1275342"
## [21] "SRR1275363" "SRR1275280" "SRR1275264" "SRR1275369" "SRR1274131"
## [26] "SRR1275351" "SRR1275256" "SRR1275355" "SRR1275258" "SRR1275260"
## [31] "SRR1275367" "SRR1275346" "SRR1274114" "SRR1275297" "SRR1275273"
## [36] "SRR1274126" "SRR1275290" "SRR1275274" "SRR1274121" "SRR1275341"
## [41] "SRR1275336" "SRR1274119" "SRR1274087" "SRR1274113" "SRR1275246"
## [46] "SRR1275267" "SRR1275289" "SRR1275360" "SRR1275358" "SRR1275352"
## [51] "SRR1274130" "SRR1275368" "SRR1275281" "SRR1275362" "SRR1275257"
## [56] "SRR1275350" "SRR1274129" "SRR1275298" "SRR1274123" "SRR1275343"
## [61] "SRR1275334" "SRR1274085" "SRR1274111" "SRR1275349" "SRR1275244"
## [66] "SRR1275344" "SRR1274088" "SRR1274116" "SRR1275243" "SRR1275301"
## [71] "SRR1275295" "SRR1275271" "SRR1274124" "SRR1275357" "SRR1275365"
## [76] "SRR1275268" "SRR1275361" "SRR1275266" "SRR1275288" "SRR1274133"
## [81] "SRR1275353" "SRR1275359" "SRR1274120" "SRR1275275" "SRR1275291"
## [86] "SRR1274086" "SRR1274112" "SRR1274118" "SRR1275340" "SRR1274115"
## [91] "SRR1275347" "SRR1274127" "SRR1275272" "SRR1275302" "SRR1275296"
## [96] "SRR1275354" "SRR1275259" "SRR1275285" "SRR1275366"
pheatmap::pheatmap(cor(exprSet))
group_list <- sample_ann$Biological_Condition
tmp <- data.frame(g = group_list)
rownames(tmp) <- colnames(exprSet)
# 组内的样本的相似性应该是要高于组间的!
pheatmap::pheatmap(cor(exprSet), annotation_col = tmp)
dim(exprSet)
## [1] 16759 99
exprSet = exprSet[apply(exprSet, 1, function(x)
sum(x > 1) > 5), ]
dim(exprSet)
## [1] 11337 99
dim(exprSet)
## [1] 11337 99
exprSet <- exprSet[names(sort(apply(exprSet, 1, mad), decreasing = T)[1:500]), ]
dim(exprSet)
## [1] 500 99
M <-cor(log2(exprSet + 1))
tmp <- data.frame(g = group_list)
rownames(tmp) <- colnames(M)
pheatmap::pheatmap(M, annotation_col = tmp)
table(sample_ann$LibraryName)
##
## GW16_1 GW16_10 GW16_12 GW16_13 GW16_14 GW16_15 GW16_16
## 1 1 2 2 2 2 2
## GW16_17 GW16_18 GW16_19 GW16_2 GW16_20 GW16_21 GW16_22
## 2 2 2 1 2 2 2
## GW16_23 GW16_24 GW16_25 GW16_4 GW16_5 GW16_6 GW21_1
## 2 2 2 2 2 1 1
## GW21_2 GW21_3 GW21_5 GW21_6 GW21_7 GW21_8 GW21+3_1
## 2 2 2 1 2 1 2
## GW21+3_10 GW21+3_11 GW21+3_12 GW21+3_13 GW21+3_14 GW21+3_16 GW21+3_2
## 2 1 1 1 1 1 2
## GW21+3_3 GW21+3_4 GW21+3_5 GW21+3_6 GW21+3_7 GW21+3_8 GW21+3_9
## 2 1 2 2 2 1 2
## NPC_1 NPC_10 NPC_11 NPC_12 NPC_13 NPC_14 NPC_15
## 2 2 2 2 2 2 2
## NPC_2 NPC_3 NPC_4 NPC_5 NPC_6 NPC_7 NPC_8
## 2 2 2 1 2 2 2
## NPC_9
## 2
可以看到,从细胞的相关性角度来看,到NPC跟另外的GW细胞群可以区分的很好,但是GW本身的3个小群体并没有那么好的区分度。
而且简单选取top的sd的基因来计算相关性,并没有很明显的改善。
但是可以看到每个细胞测了两次,所以它们的相关性要好于其它同类型的细胞。
如果计算机资源不够,这里可以先对基因进行一定程度的挑选,最简单的就是选取top的sd的基因,这里略。
dat <- dat_back
hc <- hclust(dist(t(dat)))
plot(hc,labels = FALSE)
clus <- cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝,即选择输出指定类别数的系谱聚类结果。
group_list <- as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数
## group_list
## 1 2 3 4
## 29 25 39 6
table(group_list,sample_ann$Biological_Condition)
##
## group_list GW16 GW21 GW21+3 NPC
## 1 0 0 0 29
## 2 20 3 2 0
## 3 15 8 16 0
## 4 1 0 5 0
可以看到GW16和GW21是很难区分开来的,如果是普通的层次聚类的话。
降维算法很多,详情可以去自行搜索学习,比如:
这里只介绍 PCA 和 t-SNE
dat <- dat_back
dat <- t(dat)
dat <- as.data.frame(dat)
plate <- sample_ann$Biological_Condition # 这里定义分组信息
dat <- cbind(dat, plate) #cbind根据列进行合并,即叠加所有列 #矩阵添加批次信息
dat[1:4, 1:4]
## A1BG A1BG-AS1 A2M A2M-AS1
## SRR1274090 0 0 0.000000 0
## SRR1275287 0 0 4.216768 0
## SRR1275364 0 0 0.000000 0
## SRR1275269 0 0 3.552694 0
table(dat$plate)
##
## GW16 GW21 GW21+3 NPC
## 36 11 23 29
# The variable plate (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[, -ncol(dat)], graph = FALSE)
head(dat.pca$var$coord) ## 每个主成分的基因重要性占比
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## A1BG 0.19046450 0.09601240 -0.17840553 -0.001507970 -0.0006057691
## A1BG-AS1 -0.02510451 0.29821319 0.03571804 0.020001929 -0.0105727109
## A2M 0.03403042 0.25458727 0.24264958 0.228512329 0.5414019044
## A2M-AS1 0.23140893 0.02900348 -0.07952678 0.356461354 0.1283450099
## A2ML1 -0.15776536 0.13831288 0.10065788 0.004060288 -0.0353422367
## A2MP1 -0.04068586 -0.05584736 -0.02857416 0.018287992 0.0069603680
head(dat.pca$ind$coord) ## 每个细胞的前5个主成分取值。
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SRR1274090 40.251912 -13.231641 -12.358891 -20.038100 -12.704947
## SRR1275287 1.196637 15.386256 30.566235 14.262858 -4.852418
## SRR1275364 -34.731051 -14.782146 -7.716928 7.046918 1.951473
## SRR1275269 21.760471 3.307309 17.985263 -18.382512 9.270646
## SRR1275263 -3.313968 -15.856721 8.929275 -36.358830 20.275875
## SRR1274117 59.378486 16.453551 -5.098901 56.245455 19.257598
fviz_pca_ind(
dat.pca,
#repel =T,
geom.ind = "point",
# show points only (nbut not "text")
col.ind = dat$plate,
# color by groups
#palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
# Concentration ellipses
legend.title = "Groups"
)
同样的,很明显可以看到NPC跟另外的GW细胞群可以区分的很好,但是GW本身的3个小群体并没有那么好的区分度。
因为计算量的问题,这里先选取PCA后的主成分,然进行tSNE,当然,也有其它做法,比如选取变化高的基因,显著差异基因等等。
# 选取前面PCA分析的5个主成分。
dat_matrix <- dat.pca$ind$coord
# Set a seed if you want reproducible results
set.seed(42)
library(Rtsne)
# 如果使用原始表达矩阵进行 tSNE耗时很可怕,dat_matrix = dat_back
# 出现Remove duplicates before running TSNE 则check_duplicated = FALSE
# tsne_out <- Rtsne(dat_matrix,pca=FALSE,perplexity=30,theta=0.0, check_duplicates = FALSE) # Run TSNE
tsne_out <- Rtsne(dat_matrix,perplexity=10)
plate <- sample_ann$Biological_Condition # 这里定义分组信息
plot(tsne_out$Y,col= rainbow(4)[as.numeric(as.factor(plate))], pch=19)
降维是降维,聚类是聚类,需要理解其中的区别。
降维与否,不同的降维算法选择,不同参数的选择得到的结果都不一样。
聚类也是一样,不同的算法,不同的参数。
# 前面我们的层次聚类是针对全部表达矩阵,这里我们为了节省计算量,可以使用tsne_out$Y这个结果
head(tsne_out$Y)
## [,1] [,2]
## [1,] -5.5201858 -21.8803070
## [2,] -0.2618134 -0.7269741
## [3,] 14.9423636 17.1572253
## [4,] -6.6313221 -14.8290355
## [5,] -4.6828620 -11.5136930
## [6,] -13.1192292 -18.8796564
opt_tsne=tsne_out$Y
table(kmeans(opt_tsne,centers = 4)$clust)
##
## 1 2 3 4
## 22 15 18 44
plot(opt_tsne, col=kmeans(opt_tsne,centers = 4)$clust, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
library(dbscan)
plot(opt_tsne, col=dbscan(opt_tsne,eps=3.1)$cluster, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")
table(dbscan(opt_tsne,eps=3.1)$cluster)
##
## 0 1 2 3 4
## 3 30 22 38 6
# 比较两个聚类算法区别
table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.1)$cluster)
##
## 0 1 2 3 4
## 1 0 0 0 38 6
## 2 0 16 0 0 0
## 3 1 14 0 0 0
## 4 2 0 22 0 0
library(M3Drop)
Normalized_data <- M3DropCleanData(counts,
labels = sample_ann$Biological_Condition ,
is.counts=TRUE, min_detected_genes=2000)
dim(Normalized_data$data)
## [1] 13405 97
length(Normalized_data$labels)
## [1] 97
class(Normalized_data)
## [1] "list"
str(Normalized_data)
## List of 2
## $ data : num [1:13405, 1:97] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:13405] "A1BG" "A2M" "A2ML1" "AAAS" ...
## .. ..$ : chr [1:97] "SRR1274090" "SRR1275287" "SRR1275364" "SRR1275269" ...
## $ labels: chr [1:97] "NPC" "GW21+3" "GW16" "GW21" ...
这个包设计比较简单,并没有构建S4对象,只是一个简单的list而已。
需要深入读该文章,了解其算法,这里略过,总之它对单细胞转录组的表达矩阵进行了一系列的统计检验。
fits <- M3DropDropoutModels(Normalized_data$data)
# Sum absolute residuals
data.frame(MM=fits$MMFit$SAr, Logistic=fits$LogiFit$SAr,
DoubleExpo=fits$ExpoFit$SAr)
## MM Logistic DoubleExpo
## 1 1651 1646 4033
# Sum squared residuals
data.frame(MM=fits$MMFit$SSr, Logistic=fits$LogiFit$SSr,
DoubleExpo=fits$ExpoFit$SSr)
## MM Logistic DoubleExpo
## 1 403 345 1962
DE_genes <- M3DropDifferentialExpression(Normalized_data$data,
mt_method="fdr", mt_threshold=0.01)
dim(DE_genes)
## [1] 180 3
head(DE_genes)
## Gene p.value q.value
## ABCE1 ABCE1 2.844523e-05 3.020909e-03
## ACAT2 ACAT2 7.885891e-05 6.776306e-03
## ADGRB3 ADGRB3 1.083437e-04 8.443879e-03
## AMER2 AMER2 1.917979e-05 2.255308e-03
## ANKRD44 ANKRD44 6.749811e-05 6.113596e-03
## ANP32E ANP32E 1.547604e-07 3.771932e-05
这里是针对上面的统计结果来的
par(mar=c(1,1,1,1))
heat_out <- M3DropExpressionHeatmap(DE_genes$Gene, Normalized_data$data,
cell_labels = Normalized_data$labels)
可视化了解一下找到的差异基因在不同的细胞类型的表达分布情况。
这里可以重新聚类后,针对自己找到的类别来分别找marker基因,不需要使用测试数据自带的表型信息。
cell_populations <- M3DropGetHeatmapCellClusters(heat_out, k=4)
library("ROCR")
marker_genes <- M3DropGetMarkers(Normalized_data$data, cell_populations)
table(cell_populations,Normalized_data$labels)
##
## cell_populations GW16 GW21 GW21+3 NPC
## 1 0 0 0 29
## 2 14 8 19 0
## 3 4 1 2 0
## 4 16 2 2 0
head(marker_genes[marker_genes$Group==4,],20)
## AUC Group pval
## ADGRV1 0.9707792 4 1.831217e-11
## TFAP2C 0.9451299 4 1.885637e-11
## EGR1 0.9409091 4 1.159852e-09
## PLCE1 0.9233766 4 7.676760e-11
## FOS 0.9058442 4 1.806229e-08
## SLC1A3 0.9048701 4 1.542660e-09
## AASS 0.9032468 4 3.136961e-08
## ITGB8 0.8886364 4 8.823909e-08
## BCAN 0.8844156 4 1.513264e-12
## NFIA 0.8831169 4 7.724889e-08
## PTPRZ1 0.8707792 4 3.596940e-07
## HES1 0.8590909 4 7.233415e-08
## DUSP10 0.8584416 4 7.865906e-09
## LTBP1 0.8545455 4 1.872381e-07
## LDLR 0.8487013 4 1.547451e-06
## FBXO32 0.8477273 4 3.810723e-09
## CLU 0.8467532 4 1.319982e-06
## CCND2 0.8422078 4 2.668888e-06
## JUN 0.8357143 4 8.979968e-07
## GPM6B 0.8350649 4 4.121035e-06
marker_genes[rownames(marker_genes)=="FOS",]
## AUC Group pval
## FOS 0.9058442 4 1.806229e-08
也可以针对这些 marker genes去画热图,当然,得根据AUC和P值来挑选符合要求的差异基因去绘图。
par(mar=c(1,1,1,1))
choosed_marker_genes=as.character(unlist(lapply(split(marker_genes,marker_genes$Group), function(x) (rownames(head(x,20))))))
heat_out <- M3DropExpressionHeatmap(choosed_marker_genes, Normalized_data$data, cell_labels = cell_populations)
如果遇到Error in plot.new() : figure margins too large报错,则单独将heat_out这行命令复制出来运行
通常是GO/KEGG等数据库,通常是超几何分布,GSEA,GSVA等算法。
这里就略过。
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] M3Drop_1.8.1 numDeriv_2016.8-1
## [3] dbscan_1.1-3 Rtsne_0.15
## [5] ROCR_1.0-7 gplots_3.0.1
## [7] factoextra_1.0.5 FactoMineR_1.41
## [9] cowplot_0.9.3 tidyr_0.8.1
## [11] scRNAseq_1.8.0 scater_1.10.0
## [13] ggplot2_3.1.0 SingleCellExperiment_1.4.0
## [15] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [17] BiocParallel_1.14.2 matrixStats_0.54.0
## [19] Biobase_2.42.0 GenomicRanges_1.34.0
## [21] GenomeInfoDb_1.16.0 IRanges_2.16.0
## [23] S4Vectors_0.20.1 BiocGenerics_0.28.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 RColorBrewer_1.1-2
## [3] rprojroot_1.3-2 tools_3.5.1
## [5] backports_1.1.3 R6_2.3.0
## [7] DT_0.5 KernSmooth_2.23-15
## [9] vipor_0.4.5 lazyeval_0.2.1
## [11] colorspace_1.4-0 withr_2.1.2
## [13] tidyselect_0.2.5 gridExtra_2.3
## [15] compiler_3.5.1 flashClust_1.01-2
## [17] labeling_0.3 caTools_1.17.1.1
## [19] scales_1.0.0 stringr_1.3.1
## [21] digest_0.6.18 rmarkdown_1.10
## [23] XVector_0.22.0 pkgconfig_2.0.2
## [25] htmltools_0.3.6 bbmle_1.0.20
## [27] limma_3.36.5 htmlwidgets_1.3
## [29] rlang_0.3.1 shiny_1.1.0
## [31] DelayedMatrixStats_1.2.0 bindr_0.1.1
## [33] jsonlite_1.6 crosstalk_1.0.0
## [35] gtools_3.8.1 dplyr_0.7.8
## [37] RCurl_1.95-4.11 magrittr_1.5
## [39] GenomeInfoDbData_1.1.0 leaps_3.0
## [41] Matrix_1.2-14 Rcpp_1.0.0
## [43] ggbeeswarm_0.6.0 munsell_0.5.0
## [45] viridis_0.5.1 scatterplot3d_0.3-41
## [47] stringi_1.2.4 yaml_2.2.0
## [49] edgeR_3.24.0 MASS_7.3-50
## [51] zlibbioc_1.26.0 plyr_1.8.4
## [53] grid_3.5.1 gdata_2.18.0
## [55] promises_1.0.1 ggrepel_0.8.0
## [57] crayon_1.3.4 lattice_0.20-35
## [59] locfit_1.5-9.1 knitr_1.21
## [61] pillar_1.3.1 ggpubr_0.1.8
## [63] reshape2_1.4.3 glue_1.3.0
## [65] evaluate_0.12 httpuv_1.4.5
## [67] gtable_0.2.0 purrr_0.3.0
## [69] assertthat_0.2.0 xfun_0.4
## [71] mime_0.6 xtable_1.8-3
## [73] later_0.7.5 viridisLite_0.3.0
## [75] pheatmap_1.0.10 tibble_2.0.1
## [77] beeswarm_0.2.3 bindrcpp_0.2.2
## [79] cluster_2.0.7-1 statmod_1.4.30